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This  project  involved  the  design  and  development  of  3D  regional  biogeochemical  mod¬ 
els.  The  biogeocheirdcal  models  were  designed  for  attachment  to  the  Harvard  Hierachy  of 
models,  including  the  surface  boundary  layer  (SBL)  componenl  (Walstad  and  Robinson, 
submitted,  1991)  of  the  open  ocean  baroclinic  eddy  resolving  quasigeostrophic  (QG)  model 
(Miller,  Robinson  and  Haidvogel,  J.  of  Comp.  Physics,  50:38-70,1981),  or  its  coastal  ex¬ 
tensions  (Milliff,  Math,  and  Comp,  in  Simul.,  31:541-564,1989;  and  Ozsoy,Lozano  and 
Robinson,  Ibidem,  In  Press, 1991)  (SBL/QG),  and  the  open  boundary  primitive  equation 
(PE)  model  with  hybrid  vertical  coordinates  (Spall  and  Robinson,  Math,  and  Comp,  in 
Simul.,  31:241-269,  1989).  The  selection  of  these  physical  models  is  motivated  by  the 
interest  to  study  the  time  evolution  of  biogeochemical  components  with  realistic  physiced 
fields,  as  they  all  have  data  assimilation  schemes  associated. 


The  approach  taken  in  biogeochemical  model  design  was  to  construct  first  simplified 
models  with  a  minimal  number  of  biogeochemical  components,  and  then  step-by-step  to 
introduce  further  components  and  detailed  intercomponent  interactions.  These  simplifica¬ 
tions  allow  focused  research  on  important  biological  mechanisms  in  the  presence  of  rf'alistic 
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physical  fields.  The  construction  of  regional  biogeocheinical  models  with  data  assimilation 
schemes  will  provide  a  unique  tool  to  determine  processes  and  sensitivities  operative  in  the 
real  ocean  and  to  facilitate  efficient  experimental  design  and  model  validation. 


The  lead  scientists  at  Harvard  in  the  first  stage  were  Dr.  Fred  Lipschultz  and  Dr. 
Carlos  J.  Lozano.  This  initial  effort  was  interrupted  when  Dr.  Lipschultz  took  a  position  at 
the  Bermuda  Biological  Station.  In  the  second  stage,  Dr.  Mark  Altabet  (WHOI)  undertook 
the  responsibilities  of  Dr.  Lipschultz  until  the  termination  of  this  contract.  In  the  following 
pciragraphs  we  outline  the  main  results  achieved  within  the  contract  period. 

The  simplest  biogeochemjeal  model  of  interest  is  of  the  nitrogen  cycle.  During  this 
contract  we  constructed  two  models  with  two  different  compartmentalizations  of  the  ni¬ 
trogen  components.  The  simplest  physical  model  resolving  realistically  physical  fields  of 
importance  to  biogeocheinical  processes  is  the  SBL/QG  model.  However,  in  the  model 
implementation  a  modular  approach  has  been  taken  in  order  to  facilitate  the  Insertion  of 
biogeocheinical  models  in  either  the  QG  or  PE  models  and  thus  to  allow  the  biogeochem¬ 
ical  models  to  be  driven  by  fields  produced  by  either  model.  In  the  QG  configuration 
extensive  sensitivity  studies  can  be  carried  out  for  given  sets  of  4D  physical  fields.  The 
PE  configuration  is  used  for  benchmarking  physical  processes.  Within  the  biogeochemical 
models  similarly  a  modular  approach  has  been  taken  to  allow  the  exploration  of  different 
compartmentalizations,  cycles  (e.g.,  C)  and  interaction  parameterization. 


During  the  first  stage  of  this  project  Dr.  Lipschultz  designed  and  implemented  a  3D 
biogeochemical  model  with  three  nitrogen  compartments,  namely  phytoplankton  nitrogen, 
zooplankton  nitrogen  and  dissolved  inorganic  nitrogen  that  includes  both  nitrate  and  am¬ 
monium  (PZD  model).  Growth  of  phytoplankton  is  either  light  or  nutrient  limited  using 
equations  from  Kiefer  and  Mitchell  (  Limnol.  Oceanogr.,  28:  770-776. 19S3).  Zooplankton 
grazing  are  described  by  equations  from  Franks  et  al  (  Mar.  Biol.,  91:  121-129,1986).  The 
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nutrient  field  below  the  mixed  layer  is  initialized  by  the  relationship  with  salinity  based 
on  Elrod  and  Kester  (Deep  Sea  Res.,  33:1313-1325,1986). 

At  the  time  of  his  departure  to  Bermuda,  Dr.  Lipschultz  has  completed  the  coding 
for  a  3D  version  of  the  model  attached  to  the  SBL/QG  model,  and  the  initial  stages  of 
calibration  were  underway.  Upon  arrival  to  the  Biological  Station  in  Bermuda  he  ported 
the  codes  to  a  minicomputer.  Dr.  Lipschultz  continues  to  maintain  a  close  scientific 
relationship  with  members  of  the  group;  and  his  modular  3D  implementation  facilitates 
the  construction  of  the  3D  version  of  the  biogeochemical  model  designed  and  tested  during 
the  second  stage  of  the  project. 

Dr.  Altabet’s  experience  with  nitrogen  cycles  and  his  direct  access  to  data  and  its 
interpretation  led  to  the  development  of  a  second  version  biogeochemical  model  in  terms  of 
three  nitrogen  compartments,  namely  small  nitrogen  particulates,  big  nitrogen  particulates 
and  dissolved  inorganic  nitrogen  (SBD  model).  In  this  model  big  nitrogen  particulates  sink 
due  to  gravity  providing  a  nitrogen  recirculation  path  towards  the  lower  portions  of  the 
water  column.  Dr.  Carlos  Lozano  designed,  implemented  and  carried  out  the  numerical 
calibration  of  a  one-dimensional  version  of  the  model  as  a  first  step  to  the  3D  model.  The 
essential  requirements  imposed  to  the  numerical  algorithm  were  the  ‘mass’  conservation, 
the  maintenance  of  noimegative  values  for  of  each  nitrogen  compartment  and  an  efficient 
treatment  of  the  fast  time  scale  related  to  big  particulates  sinking  rates.  Subsequently, 
Dr.  Altabet  carried  out  extensive  biogeochemical  calibration  and  sensitivity  studies  of  the 
SBD  model  in  a  parameter  range  appropriate  for  the  Bermuda  area.  Both  the  seasonal 
and  interannual  cycles  show  good  qualitative  agreement  with  observed  data,  and  the  model 
response  to  parameter  changes  is  adequate.  This  work  will  be  used  to  calibrate  the  3D 
version  of  the  model  Further  details  are  contained  in  the  appended  internal  report  which 
is  an  integral  jiart  of  this  final  report. 
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Even  though  no  immediate  publications  have  been  produced  during  the  contract  pe¬ 
riod,  substantial  progress  was  made  in  model  design  and  development  and  the  investigators 
have  accomplished  the  fundamental  basis  for  a  significant  interdisciplinary  modelling  effort 
based  on  realistic  physical  fields. 
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A  1-dimensIonal  Upper  Ocean  Biogeochemical  Model: 
Calibration  and  Simulations  for  the  Bermuda  Area 


M.  Altabet^  and  C.  J.  Lozano^ 
August  6,  1990 


1.  Introduction 

In.  this  report  we  describe  the  formulation,  numerical  implementation,  and  calibration 
of  a  1-dimensional  biogeochemical  (BGC)  model  to  be  attached  to  the  Harvard  Open  Ocean 
Surface  Boundary  Layer  Model  (Lozano,  Walstad  and  Robinson,  1990).  The  BGC  model 
has  been  constructed  with  the  following  attributes: 

l)  The  BGC  model  is  capable  of  simulating  near-surface  nitrate  fluxes,  new  produc¬ 
tion,  and  particle  flux  with  a  relatively  small  number  of  compartments  and  interactions 
(Figure  l).  A  basic  premise  for  our  effort  is  that  a  principal  influence  of  ocean  physics  on 
biology  and  chemistry  is  through  the  transport  of  nitrate  into  the  euphotic  zone.  This  flux 
combined  with  nitrogen  export  from  the  euphotic  zone  in  the  form  of  large,  fast-sinking  par¬ 
ticles  determines  population  sizes  and  particle  concentrations.  Hence  instead  of  a  classical 
NPZ  formulation  (nitrate,  phytoplankton,  zooplankton)  we  have  chosen  a  3-compartment, 
NSB  formulation  (nitrate[DIN],  small  particles  [SPN],  large  particles  [BPN])  to  reflect  the 
importance  of  the  export  term.  These  modeled  quantities  can  be  easily  compared  to  field 
data  such  as  those  found  in  Altabet  (1989,  Limnoi.  Oconnogr.).  Phytoplankton  are  an 
implicit  compohent  of  the  SPN.  Euphotic  zone  recycling  of  nitrogen  which  is  a  function 
of  zooplankton  and  bacterial  activity  is  also  not  explicit  in  this  first  version  of  the  model 
since  it  does  directly  influence  new  production. 
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2)  The  model  is  capable  of  high  vertical  resolution  in  the  upper  200  in.  In  oligotrophic 
oceans,  h’gli-light,  low  nutrient  strata  overlay  low-light,  high  nutrient  strata  giving  rise  to 
vertical  structure  in  particle  concentration  and  fluxes.  In  3-D  simulators,  nutrient  fluxes 
are  likely  to  be  sensitive  to  both  the  vertical  gradients  in  light  and  phytoplankton  and 
along  isopycnal  gradients  in  nitrate.  As  discussed  in  the  results  section,  simulation  of 
vertical  structure  is  an  important  test  of  model  validity.  High  vertical  resolution  will  also 
be  important  for  comparison  of  simulated  chlorophyll  concentration  with  remotely  sensed 
data  (and  conversely  its  assimilation  into  the  model)  since  strong  subsurface  chlorophyll 
maxima  are  persistent  features  of  the  open  ocean.  The  biogeochemical  model  (Fasham 
et  ah,  submitted  to  J.  Mar.  Res.)  which  has  been  attached  to  the  Princeton  U.  3-D 
physical  model  (basin  scale  coarse  grid)  only  considers  the  mixed  layer  without  any  vertical 
structure. 

The  uptake  of  nitrate  (DIN)  into  small  particles  (SPN)  is  the  principal  non-linear  term 
of  the  model  and  is  dependent  upon  nitrate  concentration,  small  particle  concentration, 
and  light  intensity.  The  dependencies  on  light  and  nitrate  concentration  are  hyperbolic. 
Light  levels  decrease  exponentially  with  depth  as  a  function  of  the  extinction  coefficient 
which  in  future  versions  will  be  derived  from  SPN  concentrations.  The  fraction  of  SPN 
which  is  active  phytoplankton  and  the  fraction  of  total  N  uptake  wliich  is  the  form  of 
nitrate  are  presently  fixed  parameters  but  also  will  be  variable  in  future  versions  with 
more  complex  trophic  structure.  Interconversion  between  small  and  large  (BPN)  particles 
are  simple  first-order  processes  with  fixed  parameters.  The  breakup  of  large  particles  is 
thought  to  be  as  significant  as  their  formation  in  the  ocean.  Downward  particle  flux  is 
the  product  of  BPN  concentration  and  sinking  speed  (fixed  parameter).  SPN  decay' and 
production  of  nitrate  is  considered  as  single,  first  ordered  process  inhibited  by  light  since 
nitrification  is  known  to  be  strongly  inhibited  within  the  cuphotic  zone. 

In  the  1-D  model  physical  parameters  are  defined  externally.  We  have  included  (Figure 
lb)  short  wave  radiation  I,  vertical  diffusivity  discontinuous  at  the  base  of  the  mixing  layer. 
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A  mixed  layer  of  specific  depth  and  Kj.  can  be  specified.  Below  the  mixed  layer,  is  set 
to  lower  specified  value  and  vertical  advection.  These  parameters  can  be  time  varying  in 
order  to  simulate  daily,  seasonal,  or  yearly  cycles. 

Section  2  states  the  model  governing  equations  and  some  of  it.s  integral  properties 
are  presented  in  Section  3.  The  numerical  solution  algorithms  arc  contained  in  Section  4; 
diagnosis  variables  are  introduced  in  Section  5.  Boundary  conditions  for  the  runs  described 
in  Section  6  are:  1)  no  flux  at  surface  and  2)  no  flux  at  bottom  boundary.  The  model  is, 
however,  presently  capable  of  having  fluxes  through  the  bottom  based  on  gradients  in  its 
vicinity.  Our  final  conclusions  are  drawn  in  Section  7. 
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2.  Governing  Equations 


In  the  first-order  model  Pq  we  consider  an  ocean  with  no  mean  motions  and  the 
significant  spatial  variations  occur  in  the  vertical  only.  The  ambient  effects  are  the  vertical 
light  distribution,  the  varying  time  and  space  difTusivities  —  mostly  due  to  variation  in 
mixing  layer  depth,  the  biogeochemical  fluxes  at  the  bottom,  and  vertical  p-dvection  due 
to  fluid  motions. 

We  now  establish  some  notation  and  state  an  initial  boundary-valued  problem.  Let 


n(z,t]  =  =  (SPN,  BPN,  DIN)  , 


denote  the  state  vector  dependent  upon  vertical  position  z  and  time  t.  Here  z  is  positive 
upwards  with  the  origin  at  the  water  surface. 


The  time  evolution  of  the  state  variables  is  given  by 


L(nl  =  Stii  —  M[n]n 


dn\  dn 


(z,i!]  £  D  with  -  D  —  {{z,t)  I  t  >  0,  — /i  <  z  <  0}  , 


where  the  species  local  interactions  is  described  by  the  first  term  on  the  right-hand  side 


Mfu!  = 


•^'21  0  , 


-A(z,f) 


+  ko{z,  t] 


0  0 


The  functions  A  and  ko  depend  upon  incoming  short  wave  radiation  and  here  are  given 


P+Iz 
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^ ~~  max  (  1 


I.+Id. 


[2c) 


with 


L  =  /o(0  exp[z/zo]  , 


(2ci) 


The  exchange  coefHcients  ki2,  ^21  are  constant.  The  difTusivity  tensor  is  taken  to  be 
diagonal 

0  0 


K  = 


0  0 
0  0 


(2e) 


and  only  the  second  component,  BPN,  sinks,  i.e.,  the  advection  matrix  W  is  given  by 


W  = 


1.7  0  0 

0  zo  -  rOaink  0 
0  0  lu 


^sink  ^  0  ) 


(2f) 


and  w  is  the  fluid  vertical  velocity. 


T.n  v/e  will  consider  the  initial  boundary-valued  problem  (2)  with  initial 

conditions 


n(^,0)  =  n(z);  >  0  ,  ;  =  1,2,3, 


(3) 


and  boundary  conditions 


Bofni  = 


K  ^  -  Wn 

az 


=  0 


z=0 


B/,(n]  =  K  -;^-Wn 

az 


z  =  -h 


(5) 


i.e.,  at  the  free  surface  there  are  no  fluxes  and  at  the  bottom  of  the  biologically  active 
water  column  there  are  fluxes.  Normally  we  will  assume  an  influx  of  dissolved  nitrogen 


>  0  , 


(6a) 
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into  the  domain;  no  flux  of  small  particulates 


=  0, 


(Gb) 


and  big  particulates  do  not  enter  the  system 

/(*)(<)  <0. 


The  main  objective  of  this  note  is  the  numerical  solution  of  the  initial  boundary-valued 
problem  (2),  (3]-(5)  referred  to  as  problem  Pq. 

We  introduce  the  shorthan'd  writing  n  >  0  for  >  0;  y  —  1,2,3.  A  function  V (;,.■!) 
is  an  upper  solution  to  Pq  in  D  if 


LjVj  >  0  ; 

(7  a) 

V(z,0)  >  n(z)  ; 

1 

lA 

N 

lA 

o 

(7b) 

no[v]  >  0  ; 

0  .<  <  <  imax  , 

(7c) 

n^V]  >  f  ; 

G  T  ^  fmax  • 

(7d) 

A  lower  solution  v  satisfies  (7)  with  the  inequalities  reversed. 

We  are  interested  only  on  nonnegative  solutions  (n  >  0)  of  Pq.  This  problem,  we  shall 
prove  later  on,  has  upper  and  lower  solutions;  but  in  general  the  least  upper  solution  V"  and 
the  largest  lower  solution  v*  may  not  coincide;  that  is  v*  <  V*  with  the  strict  inequality 
holding  in  some  portion  of  Domain  D.  To  ascertain  uniqueness  or  multiplicity  of  solutions 
and  the  concomiiant  solution  huG'rcaMons  it  is  necessary  to  establish  a  more  restrictive  set 
of  assumptions  on  Pq.  In  Section  4  a  numerical  algorithm  is  proposed  assuming  a  unique 
solution. 
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3.  rvjitrogeii  Content 


In  the  sequel,  the  property  of  the  matrix  M  in  (2a) 


=  0  , 

<.j  =  i 


(8) 


is  essential.  For  tlie  sake  of  clarity,  consider  first  the  special  c^rso  i.'d’)  ---  —  0,  i/d-l 

-tt'aink,  with  Ti'aink  depth  independent.  The  nitrogen  content  defined  by 


n 


0) 


j  =  i 


satisfies  the  equation 


(Oa) 


a,  /V  -- 


0 

dz 


y  y)  _  „,(=)„(: 


J  =  1 


dz 


(9b) 


In  particular  if  /v  —  K  we  can  write  the  last  equation  in  the  form 


d,N 


dz 


K 


^V1 

d'z 


0 

'Fz 


i.e.,  the  nitrogen  content  diffuses  in  the  water  column  with  difTu,sivity  K  and  has  a  sink  of 
strength  ^  (u;^"^t^‘^). 

Returning  to  (9)  we  find  that  the  total  nitrogen  content  in  the  water  column 


N{t]  =  /  A'(^,0  dz  , 

J~h 


(10a) 


satisfies  tiie  relation 


dN 

lit' 


^  ?"‘iR 

j  - 1  "  '  " 


j* 


(10b) 


•y 

I 


—  a 


The  time  evolution  of  the  total  nitrogen  content  satialies  the  orciinary  hiirerent.ia! 


equation 


=  T  /<^>(<) , 


with  initial  conditions 


rO  /  3  \ 

■v(«)  =  A„^L 


( 1 1  h) 


In  general,  --f  0  and  tiie  fluid  vertical  divergence  dio/dz  —  — div  u  is  not.  in'cessar iiy 

zero.  Indeed,  within  the  mixing  layer 


Ox  \pQfh,a]  ■  Oy  \i)ofh,u)  '  dz 


where  f  —  (r^iTy)  in  the  wind  stress,  pQ  the  density  of  the  water,  /  the  Coriolis  parameter, 
/i,„  the  depth  of  the  mixing  layer,  and  wq  is  the  interior  vertical  velocity.  Below  the  mixing 
layer  only  the  latter  term  remains.  The  total  nitrogen  content  then  is  given  by 


3  aO 


rr^-'^  dz  = 


We  shall  use  the  total  nitrogen  content  time  history  governed  by  (13)  and  (lib)  to  monitor 
the  quality  of  the  numerical  solution. 
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4.  Numerical  Model 


In  our  application  lOgink  is  two  orders  of  magnitude  larger  than  w]  that  is,  the  sinking 
process  of  big  particulates  has  a  fast  time  scale  ts  h/rysink  relative  to  diffusivity  time 
scales  4  ~  h'^/K  or  advection  time  scales  <„  ~  h/w.  In  our  first  algorithm  we  do  not 
exploit  this  time  scale  differences;  whereas  in  the  second  algorithm  we  use  abime  splitting 
technique  to  avoid  the  refined  grids  required  to  resolve  adequately  the  sinking  process. 

4.1  First  algorithm. 

(a)  Time  Step 

A  completely  implicit  time  step  algorithm  of  (2)  yields  the  two-level  scheme 

- =  0  J  >  1  ,  (14a) 

=  n',  (14b) 

where  =  n(z,  J At),  J  =  0,1,...,  and  At  is  the  time  step. 

Tlie  nonlinear  equation  (14a)  is  solved  iteratively  as  follows:  Given  initial  data 

yo  =  n-^  ,  (15a) 

and  positive  numbers  j  —  1,2,3,  compute  yj,;  p  >  1,  such  that 


y,,  ~  AtD-^+^lyp-ijy,,  -  Atd,lk-^'*'^d,ypj  +  AtWd.yp  =  +  AfN[y;,_i|y,,_i  ; 


P>  1 


where 


mu  0  0  0  mi2  0 

=  0  jn22  0  =  m2i  0  0  ,  (I5c) 

0  0  0  m^i  0  0 

and  boundary  con«-.itions 

Bo(ypl=0  B,.[y„]  =  ,  ''  (l5d) 

until 

:  i  =  1.2, 3 .  (le) 

Comment.  The  nonlinearity  in  (14a)  arises  from  the  terms  mu  and  m32;  and  the  splitting 
of  (2a) 

M  =  D  +  N 

achieves  the  uncoupling  of  species  in  the  left-hand  side  of  (14a)  allowing  fast  simple  algo¬ 
rithms. 

We  now  focus  our  attention  in  the  boundary  linear  differential  equation  problem 
(15b)-(l5d).  After  the  substitutions  we  have: 

,  (17a) 

y(2)[l  -  At(m22)-^  +  ']  -  +  +  Ai(w(^Y  +  ^d,x/^^  =  /(=)  ,  (17b) 

y(3)  _  At5.[(A'^^))-^  +  'cl.y(^)  +  Ai(u;(3))-^  +  ^c).y(^)  =  /i^)  ,  (17c) 

with 

/<»  =  (i.oy  +  ,  (18a) 
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/(2)  ^  +  Ai(mj,)^+‘>'<.Uv 


(18b) 


/f )  =  +  A<(>i>32)'^+'  .  IIS':) 

and  boundary  conditions 

(X(0):^  +  ia^(y(i))_(^„(iy  +  i5,y(i)  =  0;  at  z  =  0  (19a) 

(7v  (2y+i5^(y(2))  _  =  0  ;  at  ^  =  0  (19b) 

(7,^(3))J  +  i5^(,_^(3))_(i„(3))^+i5,y(3)  =  0;  at  z  =  0  (19c) 


(A'(i))^+i5,(y(i))  -  (u,(i))-^+i<9.y^^)  =  0  at  z  == -h  ,  (20a) 

{Ki^Y+^d.{y^^'>)  -  at  z  =  -h,  (20b) 

(A^(3))^  +  i5^(y(3))  _  +  -  -(/(3))-7  +  i  .  (20c) 


(b)  1-D  Boundary  Problem 

We  use  finite  differences  to  solve  (17)-(19).  Let  zm^i  =  -h  <  zmit-\  <  ■  ■  ■<  ^rnix  < 
zmio  =  0  be  a  partition  of  the  interval  ~h  ^  z  <  0  and  set 


dzmik  =  -zmik  +  zmh-\  i 
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Setting  3^^.  =  y{zm^f:)  we  approximate  (17)  by 


yr  I J  ) 

_  Atm  ■  ■) _ 1  +• 

^  dzmik  +  dzmik_i  dzm(.k  ^  dzm^k  +  i  ^ 


dzm^u  +  dzmtk-i  lV  ^ 


’lyy\!l,-^o]!Uy^d)--y 


0)  /^,..0)  ...(J) 


i  l^/k-1  ~  ’'^A.-fi 


=  /4'’  ="1’  +  w,i  (3'<'>)p„v  V ;  1  <  t  <  e  - 1 ,  i  =  1,2,3  ,(21) 


and  boundary  conditions 


J~  -  J'S'’!  -  5  \  =  0 ;  J  =  1, 2. 3 . 


=  J  =  1,2,3.  (23) 


This  finite  difTerence  system  (2l)-=(23)  for  each  species  can  be  cast  as  a  tridiagonal 
system  of  equations 


+  4>'>yl’l,  =  /4'>  0<k<e  ;  =  1, 2, 3 , 


aU)  ^  ^U) 
“o 


J  =  1,2,3, 


(25a) 


“  dzme,  '  2  0 


(25b) 
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=  - 


1„0). 


-W 


dzrnii  2 


0  I 


;  =  1,2,3, 


(25c) 


.U) 


2At 


dzm?.k  +  dzm^k+i 


AlL  .  „0)  ’ 

dzme,  + 


(26a) 


,0) 


2At 


Q  mz  -  ""  •" 

^  dzmtk  +  dzm^i,^i 


^+1  -xyO  ) 

^JL+l 


dzm(.k 


+  1 


J  =  1,2,3  , 


(26b) 


5O)  ^  1  _  (^0)  +  i  =  1,  2,  3  1  <  1:  <  ^  -  1  ,  (26c) 


jAJ) 

^  dzmii 


+ 1 


(27a) 


,(,■)  _  'A  1  .„(» 


“  dzmei'^  ^ 


(27b) 


The  system  of  equations  (24)  is  solved  iteratively  as  follows;  Set  e^-'^  -  d^J^  -  0,  and 


,0)  _ 


Sj) 


'  0)  .0)  +  fcO) 


(28) 


4" 


fc  =  i,...,f 


(29) 


Once  e[^\  d^P  have  been  determined  set  3^^^^  =  0,  and  proceeding  backwards  on  the  index 


k  set 


-  ^k  ^k  +  i  +  “Jt 


(30) 
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(c)  Integral  Properties 


After  multiplying  (21)  by 


5)L-  =  -  [dzm^k  +  dzmEk+\) 

u 


and  adding  from  fc  =  1  to  £  —  1  we  obtain  after  assuming  =  0, . . . ,  ^  —  1: 


where 


and 


+  ^  ^  +  AtY^  (/O))  J'+i  .  J  =  1, 2, . . .  , 

i=i 


=  Y 

j=i 


(iVO)) J  =  1  dzmCx  +  V  Skn[^^^  +  I  dzmltn^p^ . 

Jt=i 


(31) 


Equation  (31)  is  the  discrete  version  of  the  total  nitrogen  content  equation  (ll) 


(d)  Convergence  and  Positivity  of  Discrete  Solutions 
The  iterations  (28)-(30)  are  stable  if 

I'l'*  >  1  + id'’ +  6['>|  1  ;  =  1,2,3,  (32) 

and  the  nonlinear  iterations  (15)  converge  for  sufficiently  small  At  since  if  yp  = 
Cp(At)yp_i  and  ||Cij  <  1  for  At  small. 

The  numerical  calibration  consists  in  the  selection  of  the  parameters  grid  spacing 
dzmii^,  and  At  such  that  the  nonlinear  iteration  converges  with  a  reasonable  number  of 
iterations,  and  the  conservation  law  (31)  holds  within  roundoff  accuracy. 

U 


I  now  consider  conditions  for  which  the  scheme  (28)-(30)  maps  positive  functions 
onto  positive  functions.  First  it  is  convenient  to  explore  the  case  u;^^^  =  =  0,  = 

-u)3ink;  since  for  our  application  lo^ink  can  be  two  orders  of  magnitude  larger  than  the  fiuid 


advection  velocities. 


Boundary  Condition 


Fret  Surface.  The  discrete  boundary  condition  at  the  free  surface 

where  I  have  used  =  — uigink)  are  simplified  to  the  form 


=  2A'PV(’^3ink  dzm^,). 


We  ensure  that  is  nonnegative  at  the  boundary  provided 


>  1  ;  i  =  1  . 


A  sufficient  condition  for  (33)  is  to  set  the  maximum  step  size: 


dzmii  <  A'pVwink  • 


Within  the  mixing  layer  A'(')  increases  by  two  orders  of  magnitude  and  (34)  is  easy 
to  satisfy;  whereas  below  the  mixing  layer  a  condition  similar  to  (34)  requires  rather  small 

grid  sizes. 
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'  >» 


Bottom  Boundary.  At  the  bottom 


=  ~K} 


(2) 

(2)  ^<-1 


n, 


(2) 


dzm.it 


1 


(2)^(2) 


w:  'n 


1 

+  2 


simplifies  to 


n 


(2) 


2At /(2) 
^aink 


For  sufficiently  small  At,  condition  (34)  with  i  =  £  suffices  in  order  to  ensure  nonnegative 
values,  and  for  positive  flux  >  0  the  ‘adaptive’  condition 


(1  +  >  2At/^^V^3ink  , 

is  necessary  in  order  to  ensure  non  negative  n)  ■'  values.  No  such  condition  is  necessary  if 
the  flux  is  negative  at  the  boundary.  ^ 


Sufficient  conditions  for  >  0. 

From  (30)  we  can  see  that  >  0,  d  >  0,  e  >  0  implies  >  0. 

I  will  now  find  sufficient  conditions  for  the  inequalities  e  >  0,  d  >  0. 

a)  e  >  0.  If  0  <  <  1,  <  0,  <  0,  then  =  1  —  ~ 

e[^2i)  ~  At  mjj  >  1  —  —  At  mjj.  It  suffices  then  At  mjj  <  1  to  ensure  0  <  <  1. 

For  j  =  2,3  mjj  <  0;  thus  e  >  0  if: 

<  0,  4^^  <0,  At  Mu  <l  ■  '  (35) 

1  In  pur  application  we  will  often  use  the  condition 

/C)(i)  = /(3)(0)  >0  ;  t>0 

=0;  t>0. 

In  this  instance  the  adaptive  time  stepping  condition  will  not  be  necessary. 
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b)  d  >  0.  Assuming  (35)  and  >  0,  >  0  is  trivially  satisfied  for 


j  =  1,2;  whereas  for  j  =  3  it  suffices 


i  nf>  +  «r  V 


Comment  Notice  that  the  condition  (35)  for  the  coefficients  a\^\  depend  upon  the 


mesh  size  condition 


dzmli,  <  min 
j 


/('O')  /([■O') 

’  1^1+1 


whereas  the  uptake  process  limits  the  time  step  size. -For  typical  values  of  diffusivity  and 
vertical  velocities  K/iu  ~  10  meters  whereas  for  typical  sinking  speeds  /v/iOg-mk  -  i  meter. 
It  is  for  this' reason  that  we  investigate  ah  alternate  scheme  to  overcome  the  excessive  grid 
refinement  implied  by  (37). 


4.2  Second  Algorithm. 


To  overcome,  to  some  extent,  the  computational  constraints  of  the  fast  sinking  process, 
I  introduce  time  splitting  algorithms.  Setting 


C  =  W,d,n 


(38a) 


A  =  Wa^rii  —  Mn  —  5,(K5jn)  , 


(38b) 


where 


IF  =  W„  +W3; 


ty  0  0 

Wa  =  0  tu  0 

0  0  til 


0  0  0 

Wg  =  0  Ulginit  0 

0  0  0 


(38c) 
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we  write  (2)  in  the  form 


.9(11  +  Cn  +  An  =  0  . 


I  will  consider  two  algorithms.  The  first  and  simplest  algorithm  has  0[At)  conver¬ 
gence;  whereas  the  second  one  is  slightly  more  complicated,  but  has  0(At^)'  convergence. 


(a)  Time  Implicit 


,,  J+l/2  _  J 

At 


,.,.7  +  1  ,,7’  +  ! 

1  - 5 - ,._4J+i„;7+i/!  =  0 

At 


or,  in  more  detail 


pj  +  i  _  fjJ 

At 


+  =  0  , 


(39a) 


n 


T'  +  i  _ 


At 


- +  Wa^zii'^'^^  =  0  ,  (39b) 


n  =  n  , 


(39c) 


The  boundary  conditions  are 


^  0  at  ^  =  0, 


(30e) 


and  (5)  after  replacing  W  by  W„. 

In  fact  (39b)  is  (Tla)  after  we  replace  W  by  Wa  and  the  first  algorithm  in  4.1  can  be 
used  in  its  solution.  In  (39b)  the  term  Tygjnk  is  missing.  In  this  section,  I  concentrate  on 
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the  sinking  process  (39a).  Since  the  process  pertains  only  to  big  particulates  we  simplify 
(39a)  to  the  form  and 


At  "  °  ■ 


(30a') 


I  approximate  (39a’)  using  the  upstream  difference  scheme 


-(2)  „  .(2)  AhUaink  ,-(2)  ~  (2)n  (2)  ,  n  a  i 


The  solution  of  (40)  is  given  by' 


cic  ,  k  0, . . . ,  c  , 

dzmtk  ^ 


(41a) 


^0^^  ;  +  o-k)  =  k='l,...,i 


(416) 


The  scheme  (40)  has  the  property  >  0  then  >  0;  since  Ok  >  0,  „r'  >  0 
/  2\ 

implies  hj,  '  >  0;  furthermore,  (40)  is  unconditionally  stable. 
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5.  Diagnosis 


(a)  Preliminaries. 


The  ecological  system  governed  by  (2)  lies  in  the  upper  ocean  (---200  m).  The  uptake 


u (^ , t)  =  A 


/pTl^^'n^^Vimax 


is  light  limited  and  nutrient  limited.  Only  the  upper  layers  (<  100  m)  have  sufTicient 
light  during  the  day,  and  nutrients  are  essentially  supplied  from  the  lower  layers.  The 
vertical  diffusivity  has  a  strong  discontinuity  at  the  bottom  of  the  mixing  layer,  and  as  a 
consequence  time  and  space  scales  can  not  be'readily  defined  globally.  The  characterization 
of  the  behavior  of  system  (2)  via  non  dimensional  parameters  is  therefore  complex.  It  is 
advantageous  to  explore  this  ecological  system  in  reduced  subsets  of  the  parameter  space; 
in  particular  we  select  some  groupings  of  parameters  to  represent  seajonal  regimes. 

In  this  section,  I  introduce  some  non-dimensional  parameters  useful  in  exploring  sea¬ 
sonal  responses. 

The  instantaneous  euphotic  depth  is  defined  by  the  condition 


A(-/ie,0  =  Ao,  (-13) 

where  usually  logjQ(A(0,  t)/Ao)  ~  2.  The  seasonal  average  hf.  is  the  time  average  of  hf.{t) 
(over  a  period  of  4  months).  For  the  mixing  layer  depth  I  likewise  introduce  a 

seasonal  average  hn, .  It  is  to  be  noticed  that  the  daily  variability  might  be  substantial 
(e.g.,  {|hm(0  ~  ^^m|)a  ~  /irn  'H  Spring);  and  due  care  in  the  use  of  average  scalings  is 
required. 
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The  average  euphotic  depth  is  sensitive  to  water  mass  changes  and  latitude  (length  of 
day);  whereas  the  averaged  mixing  layer  depth  depends  upon  the  season  and  requires  in 
addition  to  the  mean  a  prescription  of  a  seasonal  trend,  /i',„  =  e.g.,  deepening  in  the 

Fall,  steady  in  the  Winter,  etc. 

The  seasonal  regime  can  be  characterized  therefore  by  a  triplet  h,„,  and  /i),, .  The 
time  scale  in  the  mixing  layer  is 


i  =  _ m 

•-''i  —  ) 

TU 


and  below  the  mixing  layer,  the  time  scale  is 


The  former  ~  10'*/ 10^  days  is  too  short,  whereas  the  latter  ~  10'*/ 10  days  has  a  several 
year  tune  scale;  thus  the  fraction  of  the  euphotic  zone 


=  1  - 


h,n  T  K:i'/2 


exposed  to  slow  diffusivity  characterizes  the  effect  of  mixing  iii  the  ecosystem.  Here  T  is 
the  seasonal  time  scale.  Since  the  sinking  of  large  particles  and  diffusivity  in  the  mixing 
layer  are  fast  processes;  we  seek  relevant  seasonal  trends  in  time  scales  given  by  and  the 
advcction  scale 


t„  ~  hjw  , 


('^  10^/1  days). 


Below  the  euphotic  zone  uptake  essentially  ceases  and  system  (2)  decouples  somewhat. 
There,  the  dissolved  nitrogen  (DIN)  receives  small  particulates  with  a  time  scale  = 
l/ki),  and  the  particulates  have  relaxation  time 


^  part  — 


hl2  +  ^21 
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Both  t,,,  and  /part  small  compared  to  4- 


At  a  given  depth  /i,  the  instantaneous  residence  time  is  defined  implicitly 


’-h 


n^^'‘[zj)dz  =  \  K  —X—  + 

^  07. 


or  in  terms  of  the  average 


^  J  7^^^\z,t)dz  , 


rO 


(.19) 


[A'  '—x—  + 

J  ^  riy 


dz 


z  =  -h 


(50) 


Assuming  K  ^1)—  ~  A'^  7^,  the  time  scale  for  the  residence  time  in  the  laver  above 

O  Z  T)  rti  ll  c 


=  —h  Is  given  by  where 


k  \  1 


1  1 

+ 


de  J  IJtntb 


(51) 


For  big  particulates  [j  =  2)  t„  ~  -  <<  ti  that  in  the  residence  time  is  controlled 

by  the  gravitationally  driven  sinking.  Both  small  particulate  and  uptake  depend  upon 
diffusivity  and  fluid  vertical  advcction  Deep  penetration  of  the  mixing  layer  ?/„.  small 
shortens  the  residence  time  when  the  residence  time  is  diffusion  dominated.  Notice  that 
residence  time  is  not  directly  dependent  upon  uptake  details. 

The  dimensional  analysis  above  suggest  hg  and  A  as  the  natural  space  and  time  scales 
with  derived  parameters  t/„;, 


1  -h 

MO) 


{^uis.x/  fb] 


(52) 


•  %  ^ 


"^part 


4(^12  +  ^2l) 


(53) 


Tadv  =  , 


(54) 


rU)  =  i  /  ( -L  + J-") 

^  "^adv  / 


(55) 


The  time  scales  are  known  a  priori  except  for  the  uptake  time  scale  that  depends 
upon  system  response  ((n(^^);,J. 

In  this  parameter  space  the  seasonality  enters  in  the  parameter  i]m  and  possibly  in 
the  seaisonality  of  the  upwelling/downwellir  ’’adv  =  {whelK). 
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6.  Calibration  and  Model  results 


The  model  was  initialized  with  data  typical  for  the  Bermuda  region  (Altabet,  1989). 
Parameter  sensitivities  were  tested  with  runs  of  150  to  300  days  under  steady  conditions. 
Runs  to  test  responses  to  seasonal  variations  in  mixed  layer  depth  ana  surface  light  intensity 
were  also  made.  All  our  runs  used  algorithm  4.1.  Time  step  was  0.5  hr.-  and  vertical 
resolution  was  0.5  m.  The  latter  was  found  to  be  a  necessary  condition  fc.  stability  and 
convergence  in  the  equations  governing  BPN  sinking.  The  algorithm  4.2  permits  coarser 
level  spacing  but  has  not  yet  been  implemented.  Accuracy  of  the  numerical  solutions  was 
tested  by  examining  conservation  of  total  nitrogen  over  the  run  (DIN  +  SPN  +  BPN  for 
model  domain)  which  was  achieved  in  all  cases  except  for  the  seasonal  runs.  In  these 
instances,  there  was  a  10%  loss  over  2  years  of  simulated  time  (Figure  3).  This  apparent 
loss  was  produced  by  roundoff  errors  resulting  from  the  very  high  gradients  in  BPN  within 
a  few  meters  of  the  bottom  boundary  as  a  result  of  the  no  flux  condition.  The  sensitivity 
runs  were  very  robust  with  respect  to  the  vertical  distribution  qf  properties  (Figure  2  as 
an  example).  A  strong  nitracline  was  found  between  80  and  120  m.  In  the  upper  80  m, 
DIN  concentrations  were  less  than  0.05  mol/1  and  are  of  the  order  observed  in  the  field  by 
chemiluminescent  techniques.  SPN  exhibits  a  maximum  in  the  vicinity  of  the  nitracline  and 
drops  off  sharply  below.  BPN  has  a  similar  distribution  but  its  maximum  is  20  m  deeper. 
Vertical  particle  flux  is  simply  the  product  of  BPN  concentration  and  sinking  speed  (100 
m/day  for  these  runs)  and  during  the  steady  state  portions  of  the  runs  must  equal  vertical 
DIN  flux  for  the  1-D  case.  These  general  features  and  overall  magnitude  of  quantities 
are  similar  to  observations  (Altabet,  1989).  Vertical  fluxes  are  higher  than  measured  at 
Bermuda  in  part  due  to  a  steeper  vertical  gradient  in  DIN  realized  by  the  model  and  can 
be  adjusted  with  different  initial  concentrations.  When  the  mixed  layer  depth  (mid)  is 
relatively  shallow  (40  m  in  these  cases)  the  vertical  positioning  of  the  features  was  only 
sensitive  to  the  surface  light  intensity  Iq,  the  half-saturation  constant  for  the  dependence  of 
phytoplankton  growth  on  light  (7^),  and  the  extinction  coefficient  (zq).  As  expected,  SPN 
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and  BPN  concentration  are  sensitive  to  the  first-order  constants  for  their  interconversion 
and  sinking  speed  ki2  and  /:2i-  At  low  values  for  sinking  speed  w^mk  the  ratio  of  SPNiBPN 
is  set  by  the  ratio  of  these  constants.  IncreaLsea  mixing  in  the  nitracline  produces  higher 
concentrations  of  particles  and  their  fluxes. 

Seasonal  variations  in  mixed  layer  depth  and  surface  irradiance  of  similar  phase  and 
magnitude  (modeled  for  the  present  as  a  sinusoidal  variation)  to  those  of  the  Bermuda 
region  produce  large  changes  in  concentrations  and  vertical  structure  (Figure  4).  Increas¬ 
ing  the  mixed  layer  depth  to  190  m  erases  the  subsurface  features  and  result  in  surface 
maxima  for  SPN.  Overall  levels  of  SPN  are  higher  due  to  increased  uptake  of  DIN  due  to 
its  increased  supply  to  the  euphotic  zone  by  mixing.  DIN  concentrations  are  correspond¬ 
ingly  reduced.  BPN  is  uniform  throughout  the  region  and  the  vertical  flux  of  particles  at 
100  m  have  increased  by  a  factor  of  4,  Upon  stratification,  concentrations  return  to  their 
summertime  values.  DIN  ’grows’  in  below  the  euphotic  zone  due  to  the  decay  of  particles. 
The  subsurface  maxima  intensifies  and  deepens  throughout  the  summertime  period.  The 
reproduction  of  the  temporal  pattern  for  the  2nd  year  shows  that  there  are  no  long-term 
adjustments  occurring  in  the  model.  These  features  of  the  seasonal  cycle  have  been  ob¬ 
served  at  Bermuda  except  for  the  simulated  wintertime  SPN  and  BPN  which  are  too  high 
(factor  of  2)  values  for  DIN  which  are  too  low  (factor  of  5  in  euphotic  zone).  This  in  part 
due  to  the  no-flux  condition  at  the  bottom.  In  the  real  world,  particle  flux  at  200  m  would 
draw  down  total  N  during  the  wintertime  deep  mixing  conditions  and  thus  the  amount  of 
nitrogen  available  to  support  high  particulate  concentrations  and  uptake. 
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7.  Conclusions 

This  preparatory  work  for  a  full  three  dimensional  biogeochemical  model  imbedded  in 
a  3-D  physical  model  shows  that  a  three  compartment  nitrogen  cycle  model  is  capable  of 
simulating  observed  nitrogen  distribution  in  the  upper  ocean. 

The  approach  to  our  next  stage  in  model  development  are; 

1)  Improved  light  dependency  parametrizations  Before  attachment  to  HOOM, 
the  present  BGG  model  will  be  expanded  to  include  a  dependence  on  light  from 
particle  concentrations.  Feedback  frorn  this  dependency  in  combination  with  BPN 
flux  through  the  lower  boundary  should  result  in  wintertime  simulations  closer  to 
reality.  Algorithms  to  simulate  chlorophyll  distributions  will  be  developed  that  take 

t 

into  account  photoadaptation. 

2)  Attachment  to  HOOM.  Initially,  idealized  eddies  and  jets  will  be  simulated  and 

* 

their  interaction  with  the  BGC  model  analyzed.  The  contributions  and  sensitivities 
of  specific  processes  will  be  tested  by  selectively  turning  them  off  and  on.  The  fidelity 
of  HOOM  in  simulating  isopycna!  fluxes  will  also  be  analyzed  at  this  point. 

3)  Further  1-D  model  complexity.  In  parallel  with  item  2,  more  complex  versions  of 
the  BGC  model  will  be  developed  and  tested  in  1-D  mode.  Goals  for  this  phase  of  the 
modeling  effort  arc  the  simulation  of  the  proportion  of  phytoplankton  nitrogen  to  SPN 
(now  parameterized),  recycled  nitrogen  fluxes,  and  primary  production  (in  addition  to 
new  production).  Effort  will  focus  on  adding  compartments  for  zooplankton,  bacteria, 
detritus,  and  DOM. 
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a)  Local  nitrogen  exchanges 


t**' 


b)  Physical  fields  influencing 
nitrogen  vertical  distribution 


Figure  1.  A  1-dimensional  Upper  Ocean  Biogeochemical  Model 
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Figure  2.  Example  of  a  sensitivity  run  time  history  of 
vertically  integrated  nitrogen  on  each  compartment  and 
total  one  year  scale. 
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Figure  3.  Same  as  Figure  2,  but  for  a  long  temi  run. 
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Figure  4.  Example  of  sensitivity  run  depth  versus  time 
contours  of  nitrogen  compartments  and  totals. 
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